#include <cblas.h> // 对于BLAS
#include <lapacke.h> // 对于LAPACK
#include <iostream>
#include <fstream>
#include <stdexcept>
#include <cmath> 
#include <vector>
#include "Sphere.hpp"
#include "Spline.hpp"
#include "U.hpp"

std::string directory = "src/build/";  // 目标目录

double pi = M_PI;

// 原始函数 f(x) = 1/(1 + x^2)
double F(double x)
{
    return 1 / (1 + x * x);
}

template <typename T>
double polynomial(double x, T p) {
    return p.SplineReturn(x);  
}

int main()
{
    std::vector<double> vec;
    std::vector<double> y1;

    // 第一组数据点
    for(int i = 0; i <= 9; i++) {
        vec.push_back(-4.5 + 9 * (double)i / 9);
        y1.push_back(F(-4.5 + 9 * (double)i / 9));
    }

    // 边界条件
    double fa1 = 1 / 26.0;
    double fb1 = 1 / 26.0;

    // 第二组数据点
    std::vector<double> vec1;
    std::vector<double> y2;
    for(int i = 0; i <= 10; i++) {
        vec1.push_back(-5 + 10 * (double)i / 10);
        y2.push_back(F(-5 + 10 * (double)i / 10));
    }

    // 边界条件
    double fa2 = 10 / std::pow(26, 2);
    double fb2 = -10 / std::pow(26, 2);

    // 创建 B-spline 插值器
    BSpline p2(vec1, y2, 3, "complete", fa2, fb2);

    // 插值点
    double points[7] = {-3.5, -3, -0.5, 0, 0.5, 3, 3.5};

    // 用于保存误差的容器
    std::vector<double> err1;
    std::vector<double> err2;

    // 计算第一种插值方法的误差
    std::cout << "误差 (第一种B样条):\n";
    for (int i = 0; i < 7; i++) {
        double y = USplineReturn(points[i], vec, y1, fa1, fb1);
        double error = fabs(y - F(points[i]));
        err1.push_back(error);
        std::cout << "x = " << points[i] << ", E_S(x) = " << error << std::endl;
    }

    std::cout << "\n误差 (第二种B样条):\n";
    // 计算第二种插值方法的误差
    for (int i = 0; i < 7; i++) {
        double y = polynomial<BSpline>(points[i], p2);
        double error = fabs(y - F(points[i]));
        err2.push_back(error);
        std::cout << "x = " << points[i] << ", E_S(x) = " << error << std::endl;
    }

    return 0; 
}
